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It is demonstrated numerically that smooth three degrees of freedom Hamiltonian systems which 
are arbitrarily close to three dimensional strictly dispersing billiards (Sinai billiards) have islands of 
effective stability, and hence are non-ergodic. The mechanism for creating the islands are corners of 
the billiard domain. 
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The motion of a point particle travelling with a 
constant speed inside a region D E M^, N > 2, un- 
dergoing elastic collisions at the regions 's bound- 
ary, is known as the billiard problem. Since the 
days of Boltzmann, scientists have been using var- 
ious billiard models to approximate the classical 
and semi-classical motion in systems with steep 
potentials (e.g. for studying classical molecu- 
lar dynamics, cold atom's motion in dark optical 
traps and microwave dynamics) . The invalidity of 
this approximation near certain types of trajec- 
tories is the main issue of this paper. Indeed, we 
examine this approximation in the most robust 
case of a scattering Sinai billiard (all the bound- 
ary components of the billiard are smooth, dis- 
persing, and their intersections are all oblique). 
Such billiards are known to be ergodic, hyperbolic 
and strongly mixing, thus small smooth defor- 
mations of the billiard boundaries do not change 
these properties. Nonetheless, it had been longed 
conjectured that by introducing smooth steep po- 
tentials which are close to the billiards, hyperbol- 
icity may be destroyed. In the two-dimensional 
settings, it had been proven analytically that tan- 
gent periodic orbits and certain corners produce 
stability islands for arbitrarily steep potentials, 
with precise estimates of the scaling of the islands 
size with the steepness parameter. Direct gen- 
eralization of these results to higher dimensions 
may produce non-hyperbolic behavior, but one 
would intuitively suspect that in the scattering 
case there will be always some unstable directions 
which will destroy stability. Here, we provide a 
mechanism for the creation of islands of effective 
stability (destroying both hyperbolicity and er- 
godicity) in the higher dimensional setting. We 
demonstrate numerically that the islands of sta- 
bility are created for arbitrarily steep potential 
in both two and three dimensional billiards. Fur- 
thermore, we show that the islands are created 



for an interval of steepness parameters, hence, for 
a fixed geometry, one may destroy an island by 
either making the potential steeper or softer. 



I. INTRODUCTION 

Sinai billiards are known to be ergodic and strongly 
mixinff [H m 113 . In many applications H El El 
the billiard's flow is a simplified model which imitates the 
conservative motion in a steep potential: 
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where c may be infinite. Here we always take the parti- 
cle's energy, h, to be smaller than c so that the particle 
is confined to D. An important question is whether the 
billiard and the smooth flows arc similar for sufficiently 
small e - in particular - whether the billiard's ergodicity 
property is preserved. A definite answer to such a ques- 
tion requires a well defined limiting procedure |l9l l22l| . 
For finite-range axis-symmetric potentials it was shown 
that some configurations remain ergodic H H ITtI l25l |. 
while other configurations may possess stability islands 
P, Q • Recently, it was established that in the most gen- 
eral two-dimensional settings of dispersing billiards (not 
necessarily axis-symmetric nor of finite range) the an- 
swer is definitely negative; it was proved that there are 
two mechanisms for the creation of stability islands for 
arbitrarily small e. One mechanism is a tangency - pe- 
riodic orbits or homoclinic orbits which are tangent to 
the billiard's boundary produce islands Another 
mechanism are corners - a sequence of regular reflections 
which begins and ends in a corner (termed a corner poly- 
gon) may, under some prescribed conditions, produce 
stable periodic orbits [23. In both cases it was shown 
that a two-parameter family of potentials W{q; e, a) (e is 
the softness parameter and a is responsible for a regular 
continuous change of the billiard's geometry) possesses 
a wedge in the (e, Q!)-plane, at which the Hamiltonian 
flow has an elliptic periodic orbit. This orbit limits to 
the tangent billiard orbit/ the corner polygon as e — > 0. 
Furthermore, a method for estimating the width of the 
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stability wedge in the parameter space and of the area 
of the eUiptic islands in the phase space was developed; 
for typical potentials both quantities have a power-law 
dependence on e . These findings were realized 

experimentally using cold atoms in atom-optics billiards 
In the experiments, a mixing billiard domain is 
drawn by a fast moving laser beam which encloses cold 
atoms. A small gap is opened after an initial run time, 
and the fact that the decay rate of the remaining atoms 
depends on the gap location demonstrates that the dy- 
namics is not mixing and that some of the particles are 
trapped in stability islands. The numerical simulations 
of the experiments show that islands are indeed produced 
by corner polygons 

Much less is known on the dynamics in multi- 
dimensional billiards ( iV > 3). Motivated by the Boltz- 
mann hypothesis regarding the ergodicity of hard sphere 
gas, the ergodicity property of hard-wall semi-dispersing 
billiards were extensively studied (see and ref- 

erence there in). Nowhere dispersingergodic billiards in 
with > 3 were constructed in^l^l^ls^l- In these 
papers and in [sj examples of three-dimensional semi- 
focusing billiards with mixed phase space were presented. 
Conditions under which multi-dimensional billiards with 
finite range spherically symmetric potentials are hyper- 
bolic were found in Q. A semiclassical study of three- 
dimensional Sinai billiard was presented in pol | . Recently, 
the asymptotic expansion of regular (non-tangent, away 
from corners) motion in steep multi-dimensional poten- 
tials by integrals along an auxiliary multi-dimensional 
billiard were developed |23|. In this work the geometry 
is arbitrary, and error bounds on the billiard approxima- 
tion are found. 

Here, we demonstrate numerically, for the first time, 
that islands of stability are created for arbitrarily small 
e in both two and three dimensional soft billiards. The 
ability to locate small islands of stability in the six dimen- 
sional phase space of the highly chaotic nearly-billiard 3 

d. o.f. flow may appear to be hopeless. Three technical 
innovations enable us to establish these results numeri- 
cally. The first idea is to construct a simple symmetric 
billiard, so that instead of looking for islands of stability 
in arbitrary places, we may concentrate on the proper- 
ties of a simple periodic trajectory which exists for all 
small e values by symmetry. We examine its stability 
properties by computing the monodromy matrix of the 
local return map near this orbit. Inspired by [T5ll2^ . we 
choose a trajectory which limits, as e — > 0, to the simplest 
possible corner polygon - a cord which enters a corner 
(see the bold lines in FIG^and FIG|2Jl. Furthermore, 
in the three dimensional case, by the symmetry of the 
constructed billiard, the two non-trivial pairs of eigen- 
values of the monodromy matrix are identical, and are 
thus controlled by a single parameter. The second idea 
is that by using proper rescaling it is possible to integrate 
numerically the equations of motion for arbitrarily small 

e. Indeed, if we fix the geometry and take small e val- 
ues we encounter the usual problem of stiffness near the 



boundary. On the other hand, the equivalent increase of 
the billiard domain by a similarity factor does not intro- 
duce a serious numerical problem since Viy is small in 
the domain's interior. The third idea is that the bound- 
aries of the wedges of stability in the parameter space 
may be found numerically by a continuation scheme on 
the critical eigenvalues value. Thus the stability regions 
may be found effectively and efficiently. 

II. BILLIARD GEOMETRY 

To construct concrete examples, we define the billiard 
domains as the region exterior to several spheres Ffe with 
centers at A'' and radii r'': Tk{A'',r'') = {9 G : 

N 

^{qi - Aj^)2 = {r'^f}, AT = 2 or 3. For the two di- 

i=l 

mensional case we take three circles (FIG^I. The first 
two circles (^^'^,r^'^) = (a, ±6, r) intersect at the point 
Qc ~ (rf, 0), where d[a,b,r) — a — \Jt'^ — h"^ and the 
third circle, which has a larger radius, has (^A? .r") = 
(-R - d{a,b,r),0,R) with R r > b. The angle be- 
tween the tangents to the two circles at is given by: 

0.20 = TT - acos(l - 2—), (2) 

so that when r = b these circles are tangent and a2D = 0. 
The cord 7 = {(x, y)\x £ {—d, d), y = 0} is a corner poly- 
gon: at {x, y) ~ (— d, 0) it reflects from the large circle Fa 
according to the billiard's reflection law (0i„ = (pout = 
7r/2) and at {x,y) = (d, 0) it enters a corner. We will 
study the behavior of the smooth system near this cor- 
ner polygon, thus the closing of the billiard domain away 
from this line is irrelevant here. It may be achieved by a 
union of a finite number of dispersing smooth boundaries 
which meet at non-zero angles, or by enclosing the whole 
system in a large box. For all a > the family of billiard 
tables thus defined belong to the class of Sinai billiards - 
they are mixing dynamical systems, having one crgodic 
component and a positive Lyapunov exponent for almost 
all initial conditions. 

Similarly, in the three-dimensional case, we take 
four spheres (FIG. 121314(1 . Three spheres have equal 
radii r and have equidistant centers: (A^'2,r^'2) = 
(a,5,±V36,r), {A^,r^) ^ (a,-26,0,r). These three 
spheres intersect, for r > 26, at qc = (d, 0, 0) where 
d(a, b,r) = a — \J r"^ — 46^. The fourth sphere, of radius 
i? 3> is located at a distance 2(i from the corner point: 
(A-^,r4) = (-E - d(a,6,r),0,0,i?). The angle between 
the pairs of tangent lines to the circles of intersections of 
pairs of spheres is: 

1 3 

= acos(--(l+ ^3 _ ^2/52) )) (3) 

so r = 26 corresponds to the case ol'^d = 0. Furthermore, 
the cord 7 = {(x, y, z)\x £ (— d, d), y = z = 0} is a corner 
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FIG. 1: (Color online)The billiard geometry in the 2D case. 
A cord 7 is denoted by the bold line. 




FIG. 2: (Color online)The billiard geometry in the 3D case. 
The cord 7 is denoted by the solid line. 



polygon. Here again we can close the billiard domain 
by adding a finite number of dispersing surfaces which 
intersect each other in finite angles, or by a large box, so 
that for all a > the resulting billiard domain is compact 
and dispersing. Note that if we rescale all the spheres 
and the distances between them by a fixed scale L, the 
billiards geometry will not change and the corresponding 
corner angles remain unchanged. 
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FIG. 3: (Color online)The billiard geometry in the 3D case, 
at the cross section y = 0. The cord 7 is denoted by the bold 
line. 



FIG. 4: (Color online) The billiard geometry in the 3D case: 
at the cross-section x = Xf the radius of the circles is r/ = 
\/r^ — {xf — a)^. Dotted line; Xf = a, solid line: Xf = d, 
dashed line: d < Xf < a + r. 



III. EQUATIONS OF MOTION FOR THE 
SMOOTH FLOW 



Consider the smooth motion in this region which is 
induced by the potential W{q;wo) ~ J2k=i^k{<l',wo); 
Vkiq^wo) may be taken as the Gaussian potential asso- 



ciated with the boundary component Tk 
V{Qk{q);wo) = exp(-^), where 
distance between q and the circle Tk 



Vkiq;wo) = 
Qk{q) is the 
: Qk{q) = 



\ 



N 

E( 

i=l 



[qi — A'^Y ~ f-iid "w^o is the softness parameter. 



In the cold atom experiment wq corresponds to the width 
of the laser beam [13, and V{Qk{q);wo) corresponds to 
the averaged effective Gaussian potential which bounds 
the atoms. Previously, we established that as this po- 
tential tends to a hard wall potential [wq —f 0), regular 
reflections of the smooth flow tend to those of the billiard 
|22f ? ]. By the symmetric placement of the spheres, it is 
clear that for any wq < Wq (where min^ W{q; Wq) = h), 
there exists a periodic solution ^{t,wa) — {x{t,wo),0,0) 
which limits, as wq — > to the corner polygon 7. Notice 
that studying this system for a fixed wa and a billiard 
domain which is increased proportionally by a factor L 
(so {A^,r'') {LA^ ,Lr^)), is equivalent to studying it 
in a fixed geometry with replaced by e = wq/ L. Thus, 
by increasing the domain size we may approach the limit 
e —^ without the numerical problems associated with 
the stiff limit wq 0. 



IV. NUMERICAL COMPUTATIONS 

From the analysis of we expect that the stability 
of 7(i, e, a) will depend non-trivially on both e and the 
geometrical parameter of the billiard a and that near 
Qffc = ^ islands will appear (the limit a ^ at which 



4 



■jr 















V 









■HT 




FIG. 5: (Color online) The real part of eigenvalue A at a = 
as a function of log{t) for 2D and 3D. 
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FIG. 6: (Color online) 2D. Left: real part of eigenvalue A 
(bold) at a = 0. Right: Wedges of stability in the parameter 
space. 
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FIG. 7: (Color online) 3D. Left: real part of eigenvalue A 
(bold) at a = 0. Right: Wedges of stability in the parameter 
space. See FIG|5] for phase portraits the parameter values 
corresponding to A-D. 



the billiard is not a Sinai billiard, and thus billiard orbits 
may be trapped for arbitrarily large number of reflections 
near the eorner has not been studied in jl^). We find 
that all the regions in the (a, e) plane at which islands of 
stability associated with 7(i, e,a) exist (other islands of 
stability may eo-exist), emerge from a = at some finite 
values, and converge towards (a,e) — > (f ,0). Hence, 
we first find the stability of 7(t, e, a = 0) by computing 
the eigenvalues of the monodromy matrix of the return 
map to the local cross-section at a; = for a range of e 
values. Since there is always a pair of neutral eigenvalues 
corresponding to the flow direction, for the 2d case the 
monodromy matrix has the eigenvalues {1, 1, A, j} where 
A is the largest eigenvalue which is different from 1. In the 
3d case, due to the symmetric form of the geometry, the 
spectrum is of the form {1, 1, A, -^j A, j}.(i.e. saddle-foci 
do not appear). In FIG. [S] the real part of A is shown for 
a range of e values for the 2d and 3d cases. The large os- 
cillations from positive to negative values guarantee the 
existence of intervals of e at which Re{A} G (— 1,1) - 
on these intervals A is imaginary and belongs to the unit 
circle. In the left panels of FIGEland FIGQwe present 
an enlarged segment of FIG. |S1 with a regular e scale. 
These calculations are used to find the values of e = 
at which Re{A} = ±1, where a saddle-center and a period 
doubling bifurcations occurs respectively (in the three di- 
mensional case these are double-bifurcation points due to 
the symmetry). Then, starting at (a, e) = (O, ej) , we use 
a continuation method for finding the bifurcation curves 
for a > 0, as shown in the right panels of FIGEl and 
FIGEI In the wedges enclosed by these two curves the 
periodic orbit ^{t, e, a) is elliptic, with Flouqet multipli- 
ers exp(±ja;) (in the three dimensional case each multi- 
plier has multiplicity two), and lo varies between and tt 
as the wedges are crossed. One expects that this linear 
stability will also result in nonlinear stability for most 
(non-resonant) lo values. More elaborate study of the 
resonances and the relation to the analytic predictions 
of (2^ are of interest but are beyond the scope of the 
current paper. For the two dimensional case, we verified 
that indeed the phase portraits one obtains as a wedge of 
stability is crossed are the familiar islands which appear 
near a saddle-center and a Hamiltonian period-doubling 
bifurcations (e.g. as in the Hamiltonian Hcnon map). 

In the three dimensional case, for all lo values, the 
multipliers are in 1 : 1 resonance due to the symme- 
try. For generic systems, for almost all lo values (values 
which are non-resonant with the frequency of 7(^, 6, a)), 
we expect to have non-linear stability (see e.g. [131) ■ In- 
deed, projections of the four dimensional symplectic re- 
turn map to a; = for several (a, e) values are shown in 
FIGlHl It is demonstrated that indeed inside the wedged 
region 7(t, e, a) is nonlinearly stable for the full integra- 
tion time (approximately 4000 periods). Moreover, if we 
add a sufficiently small, a-symmetric perturbation to the 
potential (e.g. V = W + 5 cos(j/ + rf) cos(z -|- jj) with 
S,r],iJ, = 0(0.0001)) we find that the effective stability 
region still persists. For the phase-space simulations we 
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FIG. 8: (Color online) 3D. Phase portraits {y,Py) at cross-section x — 0,px > for different values of a,e — 0.04, see also 
FIG|2| Notice the different scales for the first stability wedge (B-D) and the second stability wedge (A). 



use a symplectic integrator (GniCodes which keeps 

h up to accuracy of 10~^^. Thus, we can confidently de- 
tect islands with transversal kinetic energy of up to 10~^ 
(so {py,Pz) = O(10~^)). This limits our phase-space cal- 
culations to e w 0.04 - smaller values of e produce smaller 
islands and their detection via phase space plots requires 
a higher accuracy in the integration. We stress though 
that the calculations of the bifurcation curves are accu- 
rate for much smaller e values; in these calculation only 
a single return map is computed and there exists a sharp 
transition between large positive and large negative val- 
ues of the eigenvalues (see left panels of FIG. I6I7|) . so the 
existence of elliptic regimes is guaranteed. Comparing 
the 2d and 3d wedges of stability it appears that the 3d 
wedges are indeed narrower. 



V. CONCLUDING REMARKS 

While the appearance of islands in two-degrees of free- 
dom steep Hamiltonian systems is somewhat expected, 
the mechanisms for their appearance in the higher di- 
mensional settings is not as well understood (see 0, 0| 
for some generic possibilities). Furthermore, their ap- 
pearance guarantees only effective stability due to the 
possible existence of Arnold diffusion j^J. Nonetheless, 
by KAM theory, in the non-degenerate large set 

of initial conditions belongs to KAM tori and thus stay 



forever near the stable periodic orbit. Thus, the exis- 
tence of islands in the higher dimensional setting implies 
that ergodicity is destroyed independently of the possible 
leakage out of the effective stability zone after an expo- 
nentially long time. This latter possibility suggests that 
stickiness may be an interesting event also in this higher 
dimensional setting. 



Here, we propose for the first time a mechanism for the 
creation of stability islands for smooth systems which are 
arbitrarily close to strictly dispersing three dimensional 
billiards; we showed that potentials V{q; e, a) that be- 
come arbitrarily steep as e — > 0, possess wedges in the 
(e, a)-plane at which a periodic orbit is elliptic. Thus, on 
one hand, there exist one-parameter families of poten- 
tials V{q;e,a{e)) which have a stable periodic orbit for 
arbitrarily small e. Since we showed that in the wedges 
a{e) — > q;(0) > as e — > 0, it follows that these potentials 
have islands of stability even when they are arbitrarily 
close to a hard wall dispersing (Sinai) billiards. On the 
other hand, for any fixed a € (0, -f) there exists an in- 
terval of positive e values for which islands of stability 
exist. Thus, these islands may be destroyed by either 
making the potential steeper OR softer - a somewhat 
non-intuitive result. 
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